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Quantum dots in silicon are promising candidates for implementation of solid-state quantum in- 
formation processing. It is important to understand the effects of the multiple conduction band 
valleys of silicon on the properties of these devices. Here we introduce a novel, systematic effective 
mass theory of valley-orbit coupling in disordered silicon systems. This theory reveals valley-orbit 
hybridization effects that are detrimental for storing quantum information in the valley degree of 
freedom, including non- vanishing dipole matrix elements between valley states and altered interval- 
ley tunneling. 



Isolated electrons in semiconductor systems are a 
promising candidate for quantum computation because 
they exhibit excellent control and decoherence proper- 
ties pp. Much recent progress has led to demonstrations 
of both spin- and charge-based qubits in GaAs [5]-[5] and 
Si [BHH] ■ While Si has better spin decoherence properties 
than GaAs [S], silicon's nontrivial conduction band valley 
structure is a complication [10] . 

The presence of the valley degree of freedom in Si quan- 
tum dot devices can lead to difficulty in isolating a two- 
state system to use as a qubit, because valley splitting 
energies can be the same order as both Zccman split- 
tings and orbital energy spacings [T0"l - fT5] . On the other 
hand, it has been proposed to harness this valley degree 
of freedom to define noise- resistant qubits [THl EH [H]- 
Previous studies of valley states in Si have mainly fo- 
cused on an idealized picture of the valley and orbital 
physics in which the system is taken to be disorder-free, 
and hence the valley and orbital degrees of freedom are 
good quantum numbers for the system [15[ 116] . 

It has been recognized that structural disorder, such as 
atomic steps at the heterostructure interface, alloy dis- 
order, or other types of correlated randomness, can in- 
troduce new effects such as intervalley tunneling (TTl [18] • 
Furthermore, recent experimental evidence for disorder- 
influenced valley-orbit physics has been found in both 
MOS [H US and Si/SiGe systems HQ. Studying 
disorder in silicon is especially challenging, since the 
conduction band valley states couple atomic-scale dis- 
order to the micron-scale electron confinement that is 
typical of quantum dots. To analyze this problem, re- 
searchers have used computationally intensive numerical 
techniques such atomistic tight binding [55J [53J , or ana- 
lytical approaches that assume the effective mass theory 
holds with only minor corrections necessary [HI [TSj [24] . 

In this paper, we develop a systematic disorder- 
expansion technique that successfully reproduces the re- 
sults of atomistic simulations, while retaining the appeal- 
ing physical intuition and computational efficiency of ef- 
fective mass theories. Using this technique, in addition 
to the valley mixing matrix elements noticed previously, 
we identify matrix elements that correspond to valley- 
orbit hybridization, which were previously studied in an 
approximation using the two lowest energy z— states [25] . 
We also show that the presence of these matrix elements 



leads to the emergence of effects not observed in previous 
analytical treatments. In particular, we show that disor- 
der leads to finite dipole matrix elements between valley 
states, and quantitative corrections to intervalley tunnel- 
ing. Both effects are detrimental to quantum information 
stored in valley states. 

Disorder-expansion effective mass theory. — In Si quan- 
tum dots confined to a quantum well in the absence of 
disorder, one can show through tight-binding or effective 
mass theories that the energy eigenstates ^i,±(r) form 
symmetric and antisymmetric valley doublets [15] : 
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where r is the spatial position, hi is the electronic enve- 
lope function for the i th orbital, and u±k is the periodic 
part of the Bloch function located at the conduction band 
minima k = ±k z,. Here, feg = 0-82 • 27r/a is the position 
of the valley minimum, and a = 0.543 nm is the cubic 
lattice spacing in Si. 

To calculate the effects of disorder on valley states 
in Si accurately, researchers have previously relied on 
atomistic tight-binding techniques [55] [53] [55] [57] which 
are numerical and extremely expensive computationally. 
Here, we introduce a new semi-analytical technique based 
on a systematic expansion in the matrix elements of dis- 
order. This technique allows us to understand analyti- 
cally and compute accurately the effects of interface dis- 
order much faster than was previously possible. 

We consider an unperturbed problem consisting of a 
lateral, two-dimensional confinement potential V(x,y) 
that describes a quantum dot or other device, and a one- 
dimensional, vertical confinement potential U(z) that in- 
cludes the sharp interfaces, the quantum well barriers, 
and other slowly-varying components such as an applied 
electric field. Since this problem is separable, the result- 
ing wavefunction is written as ^^(r) = Fi(x,y)ipj(z), 
where i is the x-y orbital index and j is the subband 
index [TB] , 

We solve the x-y problem using the effective mass equa- 
tion, 
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where m t = 0.19 mo is the transverse effective mass in 
Si, with too the bare electron mass. The z problem is de- 
scribed similarly, with an unperturbed Hamiltonian H®. 
However, the solution method is different because in- 
cludes sharp interface potentials that couple different val- 
ley states. There are a number of well-established tech- 
niques to solve this ID problem, including augmented ef- 
fective mass treatments [151 and tight-binding techniques 
|22j . which yield solutions of the form noted in Eq. (JlJ. 
The results shown below are obtained using the ID tight- 
binding treatment. 

We now introduce disorder through the perturbation 
potential D(r). For the "bump" geometry shown in 
Fig. [I] D(r) is zero everywhere except for the small black 
region, where it has the same height as the barrier po- 
tential. We write the full Schrodinger equation, including 
the disorder potential, as 

[H° xy + H° z + D(t)] Mr) = (3) 

We solve this equation by expanding in terms of the un- 
pertubed basis set: 

oo 

Mr)= al jkFj(x,y)Mz), (4) 

where a l j k are the expansion coefficients. The problem 
then reduces to a matrix eigenvalue problem for the co- 
efficients a l j k and the energies Ei. 

The expansion described in Eqs. ([3| and Q must be 
truncated to find a numerical solution. However, the 
method is guaranteed to succeed if a sufficiently large 
basis set is used. Since matrix elements of the disorder 
potential D(r) can be large, we may need many basis 
functions to obtain quantitative accuracy. In the prob- 
lems studied below, accurate, converged solutions can be 
obtained reasonably quickly by using tens of the ipk(z) 
basis functions and tens of the Fj(x,y) basis functions, 
leading to a dense effective Hamiltonian matrix with di- 
mensions N x N, where N ^100-500. 

Application to quantum dot systems. — We now apply 
this disorder-expansion technique to calculate the eigen- 
states of a single quantum dot, and the tunneling coef- 
ficients for a double quantum dot, in the presence of a 
disordered interface. Many previous analytic treatments 
of disordered interfaces considered only the effects of val- 
ley mixing (VM) between the two low-lying valley states. 
Here, the disorder-expansion method allows us to treat 
both VM and valley-orbit hybridization (VOH), which 
describes the mixing of orbital and valley degrees of free- 
dom, and is observed only when the basis includes more 
than one orbital degree of freedom [25]. When VOH is 
significant, the electric dipole moment between the lowest 
two states, 

p = e |d 3 r(|<Mr)| 2 -|<Mr)| 2 )r, (5) 

where — e is the electron charge, can also be significant. 
In constrast, Eq. ^ yields p = when (j) and 4>i rep- 
resent pure pure orbital states, and is on the scale of the 
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Figure 1. (Color online.) Calculation of low-lying energy 
eigenstates and electric dipole moments for a 3D quantum 
dot in a quantum well with interface disorder, demonstrating 
that valley-orbit mixing induces a substantial dipole moment. 
The simulation geometry (inset) has a quantum well thickness 
of 10 nm and a barrier height of 150 meV. The parabolic quan- 
tum dot is circular, with a diameter of L = 28.3 nm. Disorder 
is introduced as a rectangularly shaped "bump" (black region) 
in the quantum well barrier, with an :r-width of W x = 2L, 
a y-width of W y = 4L, a height of a single atom, and a 
center position at (xo,yo) = (— 0.7L, — 0.7L). The top two 
curves show the components of the dipole moment p, defined 
in Eq. |Sj|, along the x (dashed) and y (dash-dotted) direc- 
tions, as a function of electric field applied along z. The 
lower set of solid curves show the lowest six energy levels in 
the quantum dot, measured relative to the ground state. At 
low fields, these form two manifolds: a lower, S-like doublet 
and an upper P-like quadruplet. Disorder introduces two dis- 
tinct effects: valley mixing (VM) between states in the same 
valley doublet, and valley-orbit hybridization (VOH) between 
states in different valley doublets. The dipole moment is typ- 
ically substantial, but it is suppressed near the VM-induced 
anticrossing indicated by the vertical dashed line. A second 
anticrossing occurs at a higher field, where the first excited 
state changes from valley-like to orbital-like, and is accompa- 
nied by a large change in the dipole moment. 

atomic lattice spacing when (j) and (f>\ are pure pure val- 
ley states. For qubit applications, a finite dipole moment 
makes the system susceptible to charge noise [2 [23 HI] • 
We consider the specific quantum dot geometry shown 
in Fig. [I] For simplicity, we choose a 2D parabolic con- 
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finement potential for the quantum dot with an energy 
level spacing of hw — 0.5 meV, corresponding to a charac- 
teristic dot size of L = ^h/{m t uj) = 28.3 nm. We choose 
the lateral dimensions of the bump perturbation to be of 
order of L, as consistent with recent stuctural character- 
ization of Si/SiGe heterostructures [30]; specifically, we 
use W x — 2L and W y = 4L. The height of the bump 
is taken to be a single atom. The quantum well width 
is 10 nm, with a barrier height of 150 meV. We com- 
pute the z-basis functions within the ID tight-binding 
method described in Ref. [52]. The full 3D calculations 
are carried out using the disorder-expansion framework 
described above, with a basis of size 5(x) x 5(y) x 30(z). 
The results are of good accuracy, as described in Ap- 
pendix [XJ 

The results of our single-dot calculations are shown in 
Fig. [l] as a function of the perpendicular electric field. 
The curves at the top of the plot p x and p y , the x and y 
components of the electric dipole moment. For the device 
specifications used here, the dipole moment is typically 
comparable to eL. Its non-monotonic dependence on per- 
pendicular electric field can be understood by examining 
the energies of the lowest six energy eigenstates of the 
quantum dot, shown in the lower portion of Fig. [T] At 
low fields, the lowest set of six energy eigenstates splits 
into two orbital manifolds. Each of these manifolds is 
further split by a small valley splitting of order 0.1 meV. 
As the field is increased, the lowest two states undergo 
successive transitions, corresponding to level anticross- 
ings: a VM transition at about 2 x 10 6 V/m, and a VOH 
transition at about 7 x 10 6 V/m. The VM anticrossing 
is caused by a competition beween two different confine- 
ment potentials: the quantum dot and the effective con- 
finement caused by disorder. The dipole moments are 
strongly suppressed at the VM anticrossing, as shown in 
Fig. [T] although the magnitude of p is never zero. We 
note that in the limit of large orbital energy spacing huj, 
the dipole moment scales approximately as l/{huj), as 
consistent with lowest order perturbation theory. There- 
fore, two methods are available to help suppress the un- 
wanted dipole moment: using smaller dots and working 
at fields corresponding to the VM anticrossing. 

We now study the impact of VOH on interdot tun- 
neling. As noted in Refs. [HI [25], structural disorder 
induces VM, so the z-component of the wavefunction is 
no longer well described by its unperturbed eigenstate. 
Since disorder varies spatially, the z-composition of the 
wavefunction will be different from dot to dot. This leads 
to intervalley tunneling, meaning that an electron can 
change valley indices when tunneling between two dots 
[17j . When quantum information is stored in the val- 
ley indices, intervalley tunneling constitutes a loss of in- 
formation. Here, we go beyond previous work |18] by 
considering tunneling between the two sides of a double 
quantum dot, and by including VOH effects using the 
disorder-expansion technique described above. 

We again adopt a simple model for interface disorder: 
a rectangular bump at the quantum well interface, as 
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Figure 2. (Color online.) The intervalley tunnel rate r y = 

(L- 1 H \R+) between two sides of a double quantum dot in the 
presence of a "bump" at the quantum well interface, as a func- 
tion of the bump position xq. Here, the height of the bump is 
one atom, L_ refers to the lowest left-localized state, and R + 
refers to the first excited right-localized state. A schematic of 
the calculation geometry is shown in the inset. In the absence 
of disorder, the ± indices refer to unperturbed valley states. 
Results are shown for bumps with widths W x = L/2, L, and 
2L in the ^-direction, and infinite widths in the y-direction. 
The solid curves are computed using the disorder-expansion 
framework described in the main text, with a basis of size 
50(x) x l(y) x 10(z), while the dashed line is an alternative 
result for W = 2L, with a basis of size l(x) x l(y) x 2(z), which 
does not admit VOH effects. For the calculations shown here, 
the dimensions of the individual dots are the same as Fig. [T] 
and we assume an electric field of 2 x 10 5 V/m applied per- 
pendicular to the quantum well. The intervalley tunneling 
rate is substantial over a wide range of bump positions and 
widths. 

shown in the inset of Fig. [2] In this calculation, we take 
the bump width to be infinite in the y direction, but 
variable in the x direction. As anticipated in Ref. |18j . 
the intra- and inter-valley tunnel rates are comparable 
when the bump width, analogous to the disorder corre- 
lation length, is comparable to the lateral widths of the 
quantum dot. We consider two quantum dots in a bi- 
quadratic potential. The individual dot potentials have 
circular symmetry, with a diameter of L = 28.3 nm, and 
an interdot separation of d = 144 nm. We define the 
states L± and R± to be the left and right-localized states, 
obtained from the lowest two eigenstates of the left and 
right individual confinement potentials. All four states 
L± and R± are computed within the disorder-expansion 
framework using a basis set of size 50(x) x l(y) x 10(z). 
(See Appendix |A] for convergence details.) 
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To calculate the tunneling, we compute the matrix 
element of the total, double dot Hamiltonian between 
left- and right-localized states. Technical details for ef- 
ficiently computing this matrix element can be found in 
Appendix [B] We run our calculation at a low applied 
electric field, F = 2 x 10 5 V/m, so that in the absence of 
disorder, the lowest two states in each dot form a valley 
doublet. This is the regime where the valley index is most 
likely a good quantum number for quantum computing. 
With no disorder, the x and z-directions are separable, 
so the intervalley tunneling term is zero: 

r_ + = (L_\H\R + )=0 . (6) 

However, the introduction of an atomic bump leads to 
significant intervalley tunneling, as shown in Fig. [2j In 
the calculation, the dot geometry is chosen such that the 
interdot, intravalley tunnel rate in the absence of disorder 
is 2 GHz, which is a typical value observed in experiments 
[5J . Over a wide range of bump positions, we confirm that 
the intervalley tunnel rate is comparable to the intraval- 
ley tunnel rate. We find that the intervalley tunnel rate 
is largest when the bump in the interface is centered over 
one of the quantum dots. The dashed line in Fig. [^cor- 
responds to only using one x-basis function, one y-basis 
function, and two z-basis functions, corresponding to the 
simple model considered in previous studies jTHj, where 
the VOH coupling is effectively turned "off." Although 
the approximate solution is qualitatively similar to the 
accurate solution, it is not quantitatively accurate. 

Our theoretical results for intervalley tunneling are in 
reasonable agreement with recent experiments in a dou- 
ble quantum dot, where tunnel rates were measured be- 
tween a (2, 1) electron occupation and two different (1, 2) 
occupations, corresponding to the ground and lowest ex- 
cited states |S]. The small energy splitting between the 
(1, 2) configurations (~ 45 fieV) is indicative a large val- 
ley component in the excited state. (Orbital excitations 
are typically larger, in the range 0.1-1 meV [51].) The 
fact that comparable tunnel rates were observed for both 
(1,2) states (2.7 GHz vs. 3.5 GHz) indicates a strong 
intervalley matrix element. 

Finally, it is interesting to compare the numerical com- 
plexity of our scheme to that of a tight binding method. 
For the double dot considered here, we achieve good ac- 
curacy with a basis of size N — 500. The computational 
bottleneck in this procedure is diagonalizing the result- 
ing N x N dense matrix, which takes a few seconds on a 
personal computer. In contrast, the number of atoms in- 
volved in a 3D tight binding calculation (excluding atoms 
outside the quantum well) corresponds to including sev- 
eral hundred million atomic sites, which requires run- 
times of many hours on modern supercomputers |27) . 

Discussion. — We have introduced a new disorder- 
expansion effective mass technique for studying disor- 
dered silicon systems. This framework provides results 
consistent with computationally intensive tight-binding 
calculations [53], while retaining the calculational sim- 
plicity and intuitive appeal of the effective mass ap- 



proach. This approach reveals additional valley-orbit 
hybridization effects, which are responsible for a non- 
vanishing dipole moment between valley states, as well 
as intervalley tunneling. 

Both valley mixing and valley-orbit hybridization are 
problematic for storing quantum information in valley 
states, since in the presence of disorder they no longer 
afford protection against charge noise, and do not have 
consistent quantum numbers between dots. We find that 
the dipole moment can be mitigated by operating the de- 
vice at a specific applied electric field, and also by making 
the dot smaller. 

The authors acknowledge useful conversations with 
A. L. Saraiva and Belita Koiller. This work was sup- 
ported in part by ARO and LPS (W911NF-12-0607), and 
by the United States Department of Defense. The views 
and conclusions contained in this document are those of 
the authors and should not be interpreted as represent- 
ing the official policies, either expressly or implied, of the 
US Government. 



Appendix A: Convergence of the disorder expansion 

In this appendix, we examine the convergence proper- 
ties of the disorder-expansion method introduced in the 
main text. To do this, in Fig. [3j we compare the per- 
formance of our disorder-expansion method in a 2D sys- 
tem to results obtained using the tight-binding method 
of Refs. [HJ |3J. Specifically, Fig. [3] (a) shows the x- 
componcnt of the electric dipole moment p, defined in 
Eq. |5]) . To be able to compare the disorder-expansion re- 
sults to tight-binding in a reasonable time, we used a 2D 
system, so we altered the disorder from the single-atom 
square described in the main text. Analogous to the sim- 
ple disorder in 3D used in the main text, here we consider 
a single-atom bump in z with x-width W x = 70.2 nm, 
corresponding to 300 atoms. This 2D tight-binding prob- 
lem can be computed in about 10 minutes on a personal 
computer and compared to the results of the disorder ex- 
pansion. As in the main text, an electric field is applied 
along the I direction with strength 2 x 10 5 V/m, the dot 
has width L = 28.3 nm, and the quantum well is 10 nm 
thick with barrier height 150 meV. 

In Fig. [3] (a), the black curve is the tight-binding result, 
while the color lines indicate different numbers of z-basis 
functions: red is 2, blue is 10, green is 20, and purple is 
40. In all cases, 5 x-basis functions were used. As more 
z-basis functions are used, the disorder-expansion results 
become more accurate (i.e., they approach the tight bind- 
ing results). In order to quantify this further, we plot in 
the inset the percent error in the left peak as a function 
of m z , the number of z-basis functions used. For large 
m z , we observe that the error falls off like TO j 1 005 ± 007 . 

Fig. [3] (b) shows the intervalley tunneling computed 
using both disorder-expansion and tight-binding tech- 
niques. As in panel (a), the black line corresponds to the 
2D tight-binding calculation, while the colored lines cor- 
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Figure 3. (a): Comparison of disorder-expansion calcula- 
tions for dipole moment with 2D tight-binding calculations. 
A bump one atom high in I with a width in x of 300 atoms 
is centered at the lateral position xq. As in the main text, 
an electric field is applied along the z direction with strength 
2 x 10 5 V/m, the dot has width L — 28.3 nm, and the quan- 
tum well is 10 nm thick with barrier height 150 meV. The 
black, solid curve is the tight- binding result, while the col- 
ored, dashed lines correspond to different numbers of z-basis 
functions: red is 2, blue is 10, green is 20, and purple is 40. 
In all cases, 5 a;-basis functions were used. Assuming the 
tight binding results reflect exact solutions, the inset shows 
the percent error in the left peak of the main plot as a func- 
tion of m z , the number of z-basis functions used. The curve 

fit for large m z indicates that the percent error scales like 
m -i.oo5±o.oo7 for large m ^ Tlie bme crQgg ^ inc j icates tne 

percent error obtained using eight z-functions from an aug- 
mented basis set described in Appendix A. (b): Comparison of 
disorder-expansion calculations for intervalley tunneling with 
2D tight-binding calculations. The disorder and system pa- 
rameters are identical to panel (a), except that two dots are 
used, and are separated by a distance d = 150.3 nm. As 
in panel (a), the black, solid line indicates tight-binding re- 
sults, while the colored, dashed lines correspond to succes- 
sively more z-basis functions. Even though tunneling is very 
sensitive to small tails of the wavefunctions along x, the 35 x- 
basis functions used here are sufficient to ensure stability such 
that the number of z-basis functions used limits the accuracy. 

responding to disorder-expansion calculations with dif- 
ferent numbers of z-basis functions used. The system 
parameters used are identical to the dipole calculation, 
except that the disorder-expansion calculations use 35 
x-basis functions in all cases, and there are two dots, 
separated by a distance d = 150.3 nm. 

As is clear from Fig. [3j the disorder-expansion tech- 



nique does converge to the tight-binding results as ex- 
pected, with the essential physics captured at a modest 
number of basis functions. However, this convergence 
can be slow, particulaly with respect to the z basis func- 
tions. This is because the perturbations we consider in- 
volve large energy scales, which effectively shift the po- 
sitions of the I energy eigenstates. We are effectively 
trying to reconstruct this shift in position by including 
many unperturbed basis states. 

We can achieve higher accuracy with fewer basis func- 
tions by tailoring our initial choice of basis to the par- 
ticular type of disorder we include in our system. For 
example, the types of disorder that we considered in this 
paper all were of the form of single-atom bumps. This 
suggests that a better z— basis would be to include eigen- 
states not only of the unperturbed quantum well, but also 
a quantum well that is narrower by the bump height. Ef- 
fectively, this means that we are supposing that the true 
z— solution will be a sum of eigenstates of the bare well 
and eigenstates of the well where the bump covers the en- 
tire system. By using this tailored basis, we can achieve 
very high accuracy with many fewer z— basis functions. 
In the inset of [3] (a) , the blue cross ( x ) indicates the per- 
cent error obtained by using this augmented basis with 
only eight z~ functions, demonstrating significantly bet- 
ter accuracy than the "brute-force" approach with sig- 
nificantly more basis functions. Therefore, a promising 
direction for future study would be to develop physics- 
informed, tailored basis sets that can help speed conver- 
gence for more general forms of disorder. 



Appendix B: Efficient computation of matrix 
elements 



The initial, unperturbed basis used for the disorder- 
expansion calculation is separable in at least z and x — y 
(and sometimes in i and y individually as well). The 
disorder perturbation mixes the eigenstates of the unper- 
turbed problem, making them no longer separable. One 
can then use these 3D states to compute matrix elements 
of desired operators directly. However, we find that 
in practice the direct computation of these 3D matrix 
elements is computationally intensive, and often takes 
longer than the disorder-expansion calculation itself. 

To bypass this bottleneck, we exploit the separabil- 
ity of our initial basis states in order to speed up calcu- 
lation of matrix elements. While this procedure is not 
strictly necessary, it enables us to speed up our calcula- 
tions greatly. We begin with the calculation of the dipole 
matrix element, Eq. ^ in the main text. For simplicity 
of presentation, we show here the computation of only 
p x ; the calculation of p y follows similarly. Recalling the 
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definition of the expansion in Eq. Q , we write 

Px = E J ' d 3 rx{a) k a] m F. ] Fi^ k -4> m ~ a jk a° lm FjFiip k il) m 

j,k,l,m 

= E " «?*«?*) J dxd v ■ xF i F i- ( B1 ) 

Here, we suppress the arguments of the F and ip func- 
tions for notational simplicity. From this, we see that to 
compute p x , we can precompute the matrix 

M !,l = J dxdy ■ xF 3 F h (B2) 

which has dimension equal to the number of x — y basis 
elements used. Then, computing p x reduces to a simple 
sum: 

P»= E { a Wi k -oP jk a° lk )M^. (B3) 

Finally, we note that in the case of an initial basis that 
is separable in both x and y, the formula simplifies even 
further, since p x is then diagonal in both y and z. 

Next, we consider the calculation of the intervalley tun- 
neling matrix elements. We rewrite Eq. ([6]) from the main 
text as 

r_+ = (L_\ (V uu (x,y) - V R (x,y)) \R+) (B4) 
+ <L_| (f + V R (x,y) + V z (z) + D(vf) \R+) , 



where Vdd is the double-dot potential in the x — y plane, 
Vr is the right-dot potential, T is the kinetic energy 
operator, and V z is the unperturbed potential along z. 
We have grouped the second term such that it forms the 
Hamiltonian corresponding to the eigenstate, which 
lets us write 

r_+ = (L_| (V uu (x,y) - V R (x,y)) \R+)+e R+ (L-\R+) , 

where e R+ is the energy eigenvalue for By using the 

same decomposition technique as we did for the dipole 
moment calculation, we can write 

T -+ = E a k a ^ 0$ + ^ + Mf ; ) , (B6) 

3,k,t 

where the matrices are defined by 

= j dxdy (Vuu(x, y) - V R (x, y)) F^F t R , (B7) 

and 

Mf d = J dxdyF^F?. (B8) 

Here, the superscript L and R denote basis functions in 
the left and right dots, respectively. As before, these 
matrices can be precomputed to increase computational 
speed, reducing the computation of r |_ to a sum. 
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